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By means of density functional theory (DFT) calculations we have quantitatively estimated the 
impact of the spontaneous polarization (SP) of the SiC(OOOl) substrate on the electronic properties 
of quasi-freestanding graphene (QFG) decoupled from the SiC by H intercalation. To correctly 
include within standard DFT slab calculations the influence of the SP, which is a bulk property, 
on a surface confined property such as the graphene’s doping, we attach a double gold layer at the 
C-terminated bottom of the slab which introduces a metal induced gap state that pins the chemical 
potential within the gap. Furthermore, expanding the interlayer distances at the bottom of the 
slab creates a local dipole moment which counters that arising from the slab’s polar character and 
allows to align the location of the graphene’s Dirac point (DP) for cubic SiC(lll) with the chemical 
potential. Thus, the DP shifts obtained for other polytypes under the same slab model become an 
almost direct measurement of the SP-induced doping. Our results confirm the recent proposal that 
the SP induces the experimentally observed p-type doping in the graphene layer which can achieve 
DP shifts of up to several hundreds of meV (or equivalently, ~ 10 13 e/cm 2 ) for specific polytypes. 

The doping is found to increase with the hexagonality of the polytype and its thickness. For the 
slab thickensses considered (6-12 SiC bilayers) an ample, almost continuous, range of doping values 
can be achieved by tuning the number of stacking defects and their location with respect to the 
surface. The slab model is next generalized by performing large scale DFT calculations where self¬ 
doping is included in the QFG via point defects (vacancy plus a H atom) thus allowing to estimate 
the interplay between both sources of p-doping (SP- versus defect-induced) which turns out to be 
essentially additive. 

PACS numbers: 73.22.Pr, 81.05.ue, 77.22.Ej 


Introduction- Epitaxial graphene (EG) on silicon car¬ 
bide (SiC) has been demonstrated to be an excellent 
material for high-performance technological applications. 
UQ On the other hand, quasi-freestanding graphene 
(QFG) on SiC obtained from EG by hydrogen interca¬ 
lation holds an even greater promise for the realization 
of graphene-based electronic devices because of its facile 
large-scale production together with excellent transport 
characteristics @-[9j] • The latter mainly arise from the ef¬ 
ficient reduction of the interaction between graphene (G) 
and SiC, otherwise strongly coupled, by the intercalated 
H layer. Although this weak interaction fully preserves 
the Dirac cones I9H12I] a p-type doping is routinely de¬ 
tected in experiments [6. 9l. 11 flM 1 61 which can be as large 
as ~ 5.5 x 10 12 cm -2 as observed by Speck et al @] for 
QFG samples on 6H-SiC(0001) or even ~ 2.0 x 10 13 cm -2 
as recently reported on 4H-SiC(0001) by Urban et al 17|, 
while for cubic 3C-SiC(lll) Coletti et al found a slight 
n-type doping 


Based on these results, among several others [7j, [8j, [If 
211, Ristein et al [22| pointed to a relationship between 
the G doping and the hexagonality of the SiC polytype 
employed as substrate, and proposed a macroscopic spon¬ 
taneous polarization (SP) doping model whereby the sub¬ 
strate’s SP creates a pseudo-charge at the surface equiva¬ 
lent to a real acceptor layer. The model explains the sign 
of the doping (since the Si-terminated surface exhibits a 


negative SP for all polytypes) as well as its dependence 
on the hexagonality of the polytype, and has found fur¬ 
ther support in several works [ill [3 [3, [3 ■ Mam- 
madov et al 2J] have very recently corroborated such 
model by means of a systematic angle-resolved photoe¬ 
mission electron spectroscopy (ARPES) study for 3C-, 
4H- and 6H-SiC QFG-systems. Two sources of doping 
are identified in this study: a surface band bending aris¬ 
ing from the bulk dopants and the SP. The former would 
be responsible for the mild n-type doping encountered in 
3C-SiC(lll) samples 3, 18, 241 while the latter would 
account for the ~ 1.5 factor between hole dopings in 
4H- and 6H-SiC samples [24]. Although not addressed 
in that work, a third main source is the self-doping in¬ 
duced by the presence of intrinsic defects in the QFG 
(vacancies and/or adatoms), thoroughly studied within 
theory [3, HH but not so well experimentally @,130. 


In this Letter we investigate the relationship between 
the SP and the graphene’s p-type doping in perfect as 
well as defected QFG for several H-passivated hexagonal 
SiC polytypes via density functional theory (DFT) based 
calculations. The importance of employing a first prin¬ 
ciples approach to this end is manifold: it represents a 
powerful predictive tool when tuning the surface density 
of e/hole carriers in G-based devices HMHH , provides a 
unique way to quantify on equal footings the competition 
between self-doping and SP and, not the least, allows to 
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explore the validity of the macroscopic dielectric theory 
when the system size is shrinked to the nanoscale, as it is 
not obvious if a direct relationship between the SP and 
the doping charge in the QFG will still hold. Neverthe¬ 
less, first principles slab calculations -aiming to model a 
semi-infinite surface- have not yet been attempted due 
to at least three non-trivial issues which need to be re¬ 
solved: (i) combining a bulk and a surface confined prop¬ 
erty (SP and G’s doping, respectively) typically requires 
rather thick slabs to achieve convergence (this is specially 
true for dielectrics with long screening lengths), (ii) an 
appropiate boundary condition at the bottom of the slab 
which pins the chemical potential, /i, within the gap re¬ 
gardless of the selected polytype is a prerequisite to make 
any differences in the doping among them meaningful 
and, (iii) the polar character of SiC(lll)/(0001) oriented 
slabs introduces an additional electric field across the di¬ 
electric which may considerably alter the final doping 
level. Additionally, the large dispersion of the tt bands 
forming the Dirac cones requires a hyperfme sampling 
of the Brillouin Zone (BZ) for an accurate estimation of 
the doping charge (see below), thus increasing consider¬ 
ably the computational time. Paradoxically, the associ¬ 
ated low G-projected density of states (PDOS) around 
the DP makes the DP shift with respect to /x, ADP, a 
highly precise gauge for the G doping, specially from the 
experimental side. 

The first goal of this work is therefore to set up a cal¬ 
culation strategy that overcomes the above fundamental 
drawbacks. After considering different slab models we 
have found that a G/H/(SiC) n /Au 2 slab with a double 
gold layer attached to the lower C-dangling bonds meets 
satisfactorily the above requirements. The model is then 
employed to calculate the QFG dopings of different poly¬ 
types as a function of the number of SiC bilayers (BLs) in 
the slab, n. Apart from the cubic 3C-SiC(lll), we have 
considered those with largest hexagonality, namely, 2H- 
, 4F1- and 6H-SiC(0001), having a stacking defect (SD) 
every two, three and four BLs, respectively. The com¬ 
bined effect of SP and defect-induced self-doping will be 
addressed at the end of this letter. 

Methods- All presented calculations have been per¬ 
formed with the GREEN code [30] and its interface to the 
DFT based SIESTA package {3l| . We employed the gen¬ 
eralized gradient approximation 32] and included semi- 


empirical van der Waals interactions 33] to account for 


the weak G/H/SiC interaction. Dipole-dipole interac¬ 
tions between spurious slab replicas were removed via 
the usual dipole-dipole corrections [sH (DDC) thus lead¬ 
ing to vanishing electric fields in the vacuum. The rest 
of the calculation parameters are described in detail in 
Ref. [3s| . The G/H/(SiC)„/Au 2 model slabs comprised a 
(2 x 2) graphene layer placed on top of n (\/3 x -\/3)i?30 o 
SiC(lll)/(0001) bilayers (BLs) with each Si dangling 
bond at the upper surface saturated by a hydrogen atom, 
while the bottom C dangling bonds were saturated by 


placing a gold plane in registry with the carbon atoms 
and a second one following an fee stacking sequence (see 
inset in Figure [ljb)). The systems were relaxed freez¬ 
ing the inner SiC BLs to bulk-like positions optimized 
independently for each polytype. In order to obtain ac¬ 
curate doping charges, large (100x100) fc-supercells were 
employed to sample the Brillouin zone [3b] . 



FIG. 1. (a) left: SV H {z) profile along a G/H/(SiC) 6 /H slab 

depicted at the center and, right: the corresponding DOS 
projected, in ascending order, on the bottom H layer, the six 
SiC bilayers, the intercalated H layer and the G; the bonds 
between C-terminated bottom of the slab and H atoms have 
been elongated by 0.45 A with respect to their optimized 
values (1.56 vs 1.11 A). Black and blue lines correspond to 
3C-SiC(lll) and 2H-SiC(0001), respectively, (b) Same as in 
(a), but a G/H/(SiC)e/Au 2 slab with gold instead of H ter¬ 
mination is shown; the C-Au and Au-Au bonds have been 
elongated by 0.55A and 0.65A respectively. The depicted 
slab geometries in each figure correspond to the 3C-SiC(lll) 
case. 


SP induced doping- Let us first address the problems 
outlined above related to the use of a slab geometry when 
trying to estimate the G doping. For the most com¬ 
mon practice of capping the lower C-terminated surface 
with another hydrogen layer these shortcomings become 
patent. In Figure Oja) we plot the 2D averaged Hartree 
potential profiles, Vh(z), and the layer resolved PDOS 
for a G/H/(SiC)e/H slab assuming a 3C- (black lines) 
or 2H-SiC (blue) stacking. Since considerable dipole mo¬ 
ments were found in both cases (AV =0.8 and 1.5 eV 
for the 3C and 2H slabs, respectively) we elongated the 
C-H bonds with respect to their optimized values (see 
Fig. [Ha)) thus reducing the potential drops by as much 
as 0.7 eV. Unfortunately, and regardless of this expan¬ 
sion, the DP remains pinned at the chemical potential 
for both slabs. The reason is the absence of any gap 
states at the bottom of the slab, so that charge neutral¬ 
ity at the graphene layer locks fj, at the DP. A natural way 
to overcome this drawback is to introduce states within 
the gap and localized at the bottom part of the slab. If 
their associated DOS is much larger than that of the G, 
then they should pin the slab’s chemical potential just as 
bulk dopants/defects fix /i in a real surface, leaving the 
G’s Dirac cones free to trap or release electrons in order 
to screen any internal fields. Attaching a double gold 
layer at the C-ended lower layer does indeed yield the 
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desired boundary condition. In Fig. QJb) we present the 
potential profiles and PDOS of such a G/H/(SiC) 6 /Au 2 
slab for 3C- and 2H-SiC. The potential drops (AP=0.3 
and 0.5 eV) were reduced by 0.2 eV after expanding the 
Au-C and Au-Au interlayer distances by large amounts 
(>0.5 A). We stress that although the expanded bottom 
geometry is not realistic, its purpose is to modify the lo¬ 
cal dipole at the bottom of the slab in order to reduce 
the total slab’s dipole moment. Inspection of the PDOS 
show that metal induced gap states (MIGSs) appear in 
both cases as a large broad peak mainly localized at the 
interface between the last SiC BL and the gold plane, 
penetrating around three BLs into the substrate. The 
MIGSs now determine the slab’s chemical potential and, 
notably, the doping in the QFG layer has almost vanished 
for the 3C case (black line), while the 2H slab (blue) dis¬ 
plays a clear p-doping with a considerable DP shift of 
around 250 meV. We have also considered attaching sin¬ 
gle and triple gold layers at the bottom of the slab to find 
that the former (latter) yields a slight n (p)-type doping 
in the G [35]. Such slabs may therefore be employed as 
model systems where bulk dopants already induce certain 
level of doping 24 1 . 




FIG. 2. (a) DOS(fc,.E) projected on G, H and first three 

SiC BLs of a G/H/SiC semi-infinite surface calculated after 
matching the Green’s function of a G/H/(SiC)i 2 /Au 2 slab 
with that of the corresponding SiC bulk. Four different SiC 
polytype surfaces are shown: 3C, 6 Hi, 4Hi and 2H. ADP is 
indicated in each plot in meV. (b) Doping of the graphene 
layer for all G/(SiC) n /Au 2 slabs considered in this work as a 
function of n, the SiC politype and the location of the SD clos¬ 
est to the surface. Left axis gives the doping surface charge 
density, 8a, and right axis the Dirac point shift, ADP=DP—/x 
(non-linear scale). Thick horizontal lines at the right indicate 
the bulk SP associated to each polytype. 


The above analysis provides clear evidence that our 
slab model allows to address the impact that the SDs 
present in the him have on the graphene layer. Fig¬ 
ure [U(a) shows a precise picture of the Dirac cones in 
the form of PDOS (k,E) for a 3C, 6 H, 4H, and 2H sub¬ 
strate. They have been calculated for a semi-infinite 
geometry after replacing the Hamiltonian of the lower 
half of the G/H/(SiC)i 2 /Au 2 slab by that of the corre¬ 
sponding bulk polytype [38[[39i]. Apart from the negli¬ 


gible value of 14 meV for cubic 3C-SiC(lll), the graphs 
clearly reveal how the DP shifts increase with the level 
of hexagonality attaining several hundreds of meV. 

We now turn our attention on how these dopings evolve 
with the him thickness. Figure[2b) illustrates the depen¬ 
dence of the graphene doping charge 6cr (left axis) and 
ADP (right axis) as a function of the slab thickness and 
the polytype considered. For the 6 H and 4H cases we 
have defined different subsets of data depending on the 
location of the SD closest to the surface, which may be 
either between the first and second BLs (4Hi and 6 Hi) 
or the second and third (4H 2 and 6 H 2 ) or the third and 
fourth ( 6 H 3 ). Except for the 3C case, which remains with 
a negligible doping for all thicknesses, the rest of plots re¬ 
veal an overall increase of Sa with n approaching their 
respective bulk SP values (indicated by thick horizontal 
lines at the right of the plot). The 6 H and 4H polytypes 
show a stair-like behavior which can be understood by 
noting that as n is increased by one, the added BL may 
or may not introduce an additional SD in the slab. In 
the former case, the slab’s SP increases and Sa shows an 
abrupt raise, whereas in the latter case, the crystalline 
region increases and hence, the depolarization is more ef¬ 
ficient leading to a slight decrease of the doping. For the 
2H case, on the other hand, Sa shows a smooth behaviour 
with no jumps due to the absence of crystalline regions. 
However, the SP for this poly type is so large that we find 
an early saturation of Sa due to the crossing of the va¬ 
lence band maximum (VBM) by the chemical potential 
at n = 12 (see Fig. [2a)), so that for larger thicknesses 
the SiC VBs will also contribute to the screening of the 
SP. 


The QFG doping dependence displayed in Fig. [2] con¬ 
stitutes one of the central results of this work as it clearly 
establishes a relationship between the SP and QFG p- 
doping also providing an estimate for the slab size re¬ 
quired to reach the bulk SP limit; at the largest thick¬ 
ness considered, n = 12 , the doping amounts to 60-80% 
of the SP while, from extrapolation, we expect that al¬ 
most a 100% should be already reached at n > 20 BLs. 
Simple electrostatic arguments dictate that the differ¬ 
ence between Sa and the SP can be entirely attributed 
to the macroscopic electric field within the dielectric, E 
(generally known as the depolarization field) via |35|: 
eo e r E = Sa + SP. 

A further remarkable issue in Fig. [2b) is the fact that 
an ample range of dopings can be obtained in an al¬ 
most continuous way by controlling the number of SDs, 
their density and their proximity to the surface G layer. 
This is in line with the wide spread of experimental val¬ 
ues reported for the graphene’s doping within the same 
polytype. For instance, for the most commonly used 
6H-SiC(0001) substrate, our calculated dopings fall in 
the 4 — 6 e x 10 12 /cm 2 range, in excellent agreement 
with a number of experimental results varying from 2.0 
to 6.2 xlO 12 cm -2 [ 6 |, [jjlj 21, 241. Smaller Sa values 
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FIG. 3. Top (a) and side (b) view of the geometry corre¬ 
sponding to a vacancy embedded in an (8 x 8) supercell of 
G/H/2H-SiC(0001)/Au2. In (a) only graphene atoms and the 
H atom saturating the dangling bond are shown. DOS(fc,i?) 
projected on the first layers of the surface for 3C-SiC(lll) (c) 
and 2H-SiC(0001) (d). Doping levels are indicated in meV. 


(~ 2 ex 10 12 /cm 2 ) have also been reported [1, Q 
which, in light of our results, might be associated with 
defects or impurities/dopants either at the G layer or 
in the bulk [241 . or even to a cubic termination at the 
surface 14], On the other hand, few values of 5a have 
been reported for a 4H-SiC(0001) substrate: Mammadov 
et al 24] obtained 6.9 and 8.6 e x 10 12 /cm 2 for an n- 
type doped and a semi-insulating substrate, respectively, 
which is in very good accordance with our values for the 
4 H!/2 at n > 12, while the Hall measurements of Ur¬ 
ban et al 171 yielded charges of 15 — 20 e x 10 12 /cm 2 
which signifcantly exceed the bulk SP, again suggesting 
the presence of other sources of doping. Unfortunately, 
we are not aware of any QFG experiments carried out on 
a 2H surface. 


Competition between SP and self-doping- Finally, and 
in order to reach a general picture of the doping in QFG 
systems, we incorporate defects in the G layer within our 
gold terminated slab model. Among the various types of 
point defects studied in the literature [l2], [25|,|4C)| we con¬ 
sider the Jahn-Teller distorted vacancy structure where 
two C dangling bonds establish a bond among them in¬ 
ducing a pentagonal structure, while the third C dan¬ 
gling bond is saturated by a H atom (see Figure EUa)). 
We make this choice because this defect structure hardly 
alters the Dirac cones and, at the same time, yields a 
considerable p-type doping 12], Due to computational 
limitations, the defect is embedded in an (8x8) G super¬ 
cell corresponding to a somewhat large concentration of 
0.8% (typical experimental values obtained for EG sam¬ 
ples are of the order of 0.006%), while the number of SiC 
BLs was set to n=6 , leading to a total of 848 atoms in 
the unit cell (see Fig. [U]b)). 


PDOS (k,E) maps are shown in Figs. [3]]c) and (d), 
for the 3C and 2H polytypes, respectively. Although 
a band gap of ~150 meV appears between the Dirac 
bands in both cases due to the large defect concentra¬ 
tion, it is still possible to estimate their doping levels. 


For the 3C cubic case, a considerable hole concentration 
of 20.2 e x 10 12 /cm 2 (ADP=208 meV) is induced in the 
7 r bands in order to compensate the electron charge ac¬ 
cumulated in the vicinity of the defect [l2j]. As expected, 
the doping is even larger for the 2H slab, ADP=300 meV 
or 31.3 e x 10 12 /cm 2 , since the SP further contributes to 
the doping. Most importantly, the difference in 5a be¬ 
tween both polytypes is 11 e x 10 12 /cm 2 , which coincides 
very nicely with the SP-induced doping of 10 ex 10 12 /cm 2 
found for undefected 2H slabs at n = 6 (blue line in 
Fig- IHb)) indicating that the SP- and self-doping mech¬ 
anisms are basically additive. This is not a trivial result 
as it is not clear a priori if the self-doping mechanism 
will affect the depolarization field E which, as mentioned 
above, plays an important role in the final G doping. Our 
results thus suggest that the self-doping, already very 
large on its own, hardly alters the incomplete compensa¬ 
tion of the SP pseudo-charge by the G’s n bands. 

In summary, we have presented a methodology within 
the framework of standard DFT slab calculations that 
accounts for the contribution of the SiC substrate’s SP 
to the QFG doping. The scheme relies on the pinning of 
the slab’s chemical potential at the bottom of the slab via 
the creation of MIGSs states after saturating the lower 
C-dangling bonds with a gold bilayer. One may addition¬ 
ally incorporate self-doping contributions due to defects 
or adsorbates in the G layer. Application to QFG on 
6 H-, 4H- and 2H-SiC(0001) substrates indicates that full 
compensation of the SP by the G doping should occur at 
thicknesses of n > 20 BLs, while thinner slabs yield an 
ample range of dopings depending on the polytype and 
the precise termination of the surface. Other sources 
of doping, such as self-doping in the G layer, or bulk 
dopants, may also be incorporated into the calculations. 
Apart from the obvious applicability of the analysis to 
ultrathin SiC films, the scheme should also work satis¬ 
factorily in other dielectrics exhibiting an SP although 
the nature of the metallic layer and the interlayer ex¬ 
pansions required to minimize the surface dipole will in 
general need to be tuned for the specific system. 
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